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The  Equation  of  State  for  Ammonia* 
by 

Lester  Haar  and  John  Gallagher 

1.  Introduction 

This  report  has  been  prepared  in  answer  to  several  urgent  requests  for 
detailed  computer  programs  for  the  thermodynamic  properties  of  ammonia.  In 
it  we  present  a description  of  an  extensive  correlation  of  the  thermodynamic 
properties  of  ammonia  published  elsewhere'*’,  along  with  basic  equations  for 
the  description  of  the  properties  of  ammonia  and  detailed  listings  of 
computer  programs  based  on  these  equations.  With  this  information, 
thermodynamic  properties  can  be  calculated  for  temperatures  ranging  from 
the  triple  point  temperature  to  about  5/3  the  critical  temperature  and  for 
pressures  ranging  from  the  dilute  gas  to  about  8000  bars.  The  reference 
state  for  all  properties  is  the  ideal  gas  at  zero  kelvin.  The  physical 
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constants  used  are  consistent  with  those  recommended  by  Cohen  and  Taylor  . 

3 

The  mass  of  a mole  of  ammonia  was  taken  to  be  17.0306  grams  . 

2.  The  derivation  of  the  thermodynamic  surface. 

In  the  correlation,  selected  P,p,T  data  were  fitted  to  an  analytic 
equation  of  state  in  the  least  squares  sense.  A detailed  discussion  of  the 
process  of  data  selection  and  of  selection  of  parts  of  the  analytic  equation 
is  contained  in  reference  1.  The  equation  so  obtained  contains  the  pressure 
as  a 44  term  double  power  series  function  of  temperature  and  density. 

The  resulting  equation  can  be  used  to 

* This  work  was  supported  by  the  Office  of  Standard  Reference  Data  of  the 
National  Bureau  of  Standards. 
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reproduce  all  the  available  P,p,T  experimental  data  as  well  as  to  produce 
limited  extrapolations  (based  on  thermodynamic  arguments)  of  the  surface 
into  important  regions  where  data  are  sparse.  The  range  of  the  equation 
is  bounded  at  low  temperatures  by  the  triple  point  temperature  (195. 48K) 
and  the  melting  curve  for  the  liquid,  and  at  high  temperatures  by  the 
isotherm  at  750K  (which  is  approximately  5/3  the  critical  temperature) . 
The  pressure  range  extends  to  8,000  bar.  No  mathematical  constraints 
were  imposed  on  the  equation,  and  only  P ,p ,T  data  were  used  in  the 
least  squares  fit  with  the  sole  exception  that  values  of  vapor  pressure 
were  explicitly  employed  in  regard  to  satisfying  conditions  of  phase 
equilibrium. 
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Following  Keenan  et  al.  the  Helmholtz  free  energy  function  was 
represented  in  reference  1 as  the  sum  of  two  terms,  the  first  being  the 
contribution  from  the  equation  of  state  and  the  second  a function  of 
temperature  only  referring  to  the  properties  of  the  ideal  gas.  Thus  the 
Helmholtz  free  energy  was  expressed 


A(p , T)  = A(p , T)  + A°(T), 


(1) 


where  A°(T)  is  the  contribution  of  the  ideal  gas.  T is  the  absolute 
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temperature  in  kelvins  and  p the  density  in  grams  per  cm  . A quantity 


Q(p,T)  was  defined  by 


A(p,T)  = RT  [£n  p + pQ  (p,T)] 


(2) 


2 

Since  P = p 9A/dp, 


Eqs.  (1)  and  (2)  yield 

P = pRT  [1  + pQ  + p2  9Q/ 9p ] . 


(3) 
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Note  that 


Q (p=0)  = B2, 


(3a) 


where  is  the  second  virial  coefficient.  The  form  chosen  for  Q 


was 
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(4) 


i-1  j-1 


500 


where  t = — — , x^  = 1.2333498  and  R is  the  gas  constant.  Eq.  (3)  for 

the  pressure  was  fitted  in  the  least  squares  sense  to  the  experimental 

P,p,T  data.  The  results  of  this  fit  are  values  for  the  constants  a.. 

ij 

listed  in  table  1.  (a_  = ® f°r  -*->3  pairs  not  listed  in  table  1)  . 

By  differentiation  of  eq  (1) , the  various  thermodynamic  quantities  can 
be  expressed  in  terms  of  Q: 
the  entropy, 


S(p,T)  = R[£n  p + pQ  - px  3Q/3x]  + S°(T), 


(5) 


the  internal  energy, 


E(p,T)  = RpTx3Q/3x  + E°(T), 


(6) 


the  constant  volume  heat  capacity, 


Cv(p,T)  = - Rpx2  32Q/3x2  + Cv°, 


(7) 
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-6.453022304053 

-13.719926770503 

-8.100620315713 

-4.880096421085 

-12.028775626818 

6.806345929616 

8.080094367688 

14.356920005615 

-45.052976699428 

-166.188998570498 

37.^08950229818 

-40.730208333732 

1.032994880724 

55.843955809332 

492.016650817652 

1737.835999472605 

-30.874915263766 

71.483530416272 

-8.948264632008 

-169.777744139056 

-1236.532371671939 

-7812.161168316763 

1.779548269140 

-38.974610958503 

-66.922858820152 

-1.753943775320 

208.553371335493 

21348.946614397509 

247.341745995422 

299.983915547501 

4509.080578789798 

-37980.849881791548 

-306.557885430971 

24.116551098552 

-9323.356799989199 

42724.098530588371 

161.791003337459 

-507.478070464266 

8139.470397409345 

-27458.710626558130 

-27.821688793683 

298.812917313344 

-2772.597352058112 

7668.928677924520 
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the  enthalpy  function, 


H(p,T)  - RT[pQ  + p2  3Q/3p  + pT  3Q/3T  + 1]  + E°  (T)  , (8) 


the  constant  pressure  heat  capacity, 


Cp  (p,T) 


C + R* 
v 


a 

6 


t 


(8a) 


where. 


a - 1 + pQ  + p2  3Q/3p  - px  3Q/3t  - p2i  32Q/3x3p  , 


and 


8-1+2  pQ  +4  p2  3Q/3p  + p32Q/3p2, 
the  heat  capacity  for  the  saturated  fluid  , 

T OP/9T)  dPc 

s ^ OP/ap)T  dT  > (9) 

where  Pg  is  the  vapor  pressure  of  the  liquid  • 

S°,  E°  and  Cv°  are  the  corresponding  contributions  obtained  from  A°(T): 
A°(T)  - (G°  - E° ) - RT(1  + In  4.8180  T)  , 
s°  - - d?  A°(I)- 

Eo.Ao(I)  .TdA°(T),  (10) 

ai 

Cy°  - - T d2  A°(T)/dT2 » 
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where  G°  is  an  analytic  representation  of  the  results  reported  by  Haar5  for 
the  properties  of  the  ideal  gas  state, 


G - E 

o 

RT 
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= a^nTf^  a.  T1"3, 
i=2 


(10a) 


O 

where  Eq  is  the  energy  for  the  ideal  gas  at  0 K.  The  coefficients  a^  are 
listed  in  table  2. 


Table  2 - 


al 

m 

-3.872727 

a2 

m 

.64463724 

a3 

m 

3.2238759 

*4 

m 

-.0021376925 

a5 

m 

.86890833xl0"5 

a6  • 

m 

-. 24085 149xl0“7 

a7 

m 

.36893175xl0”10 

a8 

a9 

m 

r*.  35034664xl0~13 

m 

.2056303xl0”16 

aio 

- 

— .6853420xl0-20 

all 

- 

.9939243xl0"24 

The  heat  capacity  values  and  the  other  thermodynamic  functions  calculated 
from  Eq. (10a)  for  the  temperature  range  100  K < T < 1000  K agree  with  those 
tabulated  in  reference  5 to  within  the  accuracy  of  those  values. 
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Though  Eqs.  (l-10a)  are  complete,  It  is  necessary  to  introduce  the 

Gibbs  phase  conditions  in  order  to  calculate  the  properties  for  the 

coexisting  phases.  However,  it  was  shown  by  Haar  and  Gallagher'*’  that 

almost  negligible  error  results  if  an  explicit  relation  is  used  for  the 

vapor  pressure,  P , obtained  by  separately  fitting  the  saturated  vapor 

pressure  data  of  Cragoe  for  the  range  from  the  triple  point  to  373K  and 
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of  the  mean  of  the  data  of  Beattie  and  Lawrence  and  of  Keyes  above  373K. 
The  equation  so  obtained  is 

loge  = Y"  (AX  + Bx3/2  + Cx5/2  + Dx5) 
c 

where  X = 1 - T/T 

c 

P = 111.85  atm. 
c 

T = 405. 4K 
c 

A = -7.296510 
B = 1.618053 
C = -1.956546 
D = -2.114118 

Eq.  (10b)  and  Eq.  (3)  define  the  coexisting  phases  for  this  report 

and  all  properties  over  the  temperature-pressure  range  of  the  thermodynamic 

surface,  subject  to  the  low  temperature  boundary  of  the  melting  solid. 

The  thermodynamic  surface  is  consistent  with  the  following  values  for 
the  parameters  at  the  triple  point, 

T = 195. 48K 
P = .06063  bar 

pg  = .00006382  g/cm3,  = .73374  g/cm3, 

and  at  the  critical  point, 

T = 405. 4K 
c 

P 113.04  bar 
c 

P = .2350  o/cm 


8 


The  relationship  between  the  pressure  and  temperature  of  the  melting 

solid  was  calculated  by  means  of  the  Clapeyron  equation.  For  the  latent 
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heat,  the  value  reported  by  Overstreet  and  Giauque  was  used;  for  the 
specific  volume  for  the  solid  at  the  normal  melting  point,  the  value 
reported  by  McKelvey  and  Taylor'*'*’*  was  used,  and  for  the  corresponding 
specific  volume  of  the  liquid  the  value  reported  by  Cragoe  and  Harper'*’’*' 
was  used. 

The  Clapeyron  equation  is  the  relation 


dT 

T 


(ID 


where  the  quantities  u’  and  u are  the  specific  volumes  of  the  liquid 
and  solid,  respectively,  and  L is  the  latent  heat  of  fusion.  From  the 
above  data,  the  quantity  — — - — — - 4x10  atm  and  Eq.  (11)  can  be 

i-< 

integrated  to  yield 


T = Tg  exp [4x10  5(P  - P )], 


(12) 


where  Tg  and  Pg  are  the  triple  point  values.  (The  differences  between  the 
triple  point  values  and  those  of  the  normal  melting  point  are  negligible.) 
Also,  since  Pg  at  the  triple  point  is  very  small,  the  relationship  can  be 
simplified  to 


T = 195.48  exp  {4x10  ^ P (atm)}. 


(13) 
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3.  Computer  programs. 

The  properties  calculated  from  the  equations  presented  in  this  report 
have  been  compared  with  the  various  thermodynamic  measurements  reported  for 
ammonia  in  the  comprehensive  review  in  reference  1.  In  almost  all  cases 
the  agreement  is  within  the  experimental  accuracy  of  the  data.  Comparisons 
of  results  with  P,P,T  data  are  given  here  in  figures  2 and  3.  It  was 
established  in  reference  1 that  for  most  of  the  vapor  phase  and  for  the 
coexisting  phases,  the  calculated  enthalpies  are  accurate  to  within  0.1%. 

To  facilitate  the  application  of  these  results  we  present  in 
Appendices  I through  IV  computer  programs  in  FORTRAN  with  which  the  various 
properties  of  ammonia  can  be  calculated.  Four  such  programs  are  included: 

The  computer  program  in  Appendix  I is  a general  program  for  all  the 
thermodynamic  properties  discussed  in  section  2 of  this  report,  including 
the  program  for  the  calculation  of  Q of  Eq.  (4)  and  its  derivatives  with 
respect  to  temperature  and  density  which  are  used  by  most  of  the  other 
programs.  The  independent  variables  are  either  pressure  and  temperature  or 
density  and  temperature. 

The  computer  programs  in  Appendix  II  refer  only  to  certain  properties 
for  the  coexisting  phases,  including  the  latent  heat  of  vaporization,  the 
vapor  pressure  of  the  saturated  liquid,  the  densities  of  the  saturated 
vapor  and  of  the  saturated  liquid  and  the  heat  capacity  of  the  saturated 
fluid.  The  dependent  variable  for  these  is  either  the  temperature  or  the 
corresponding  pressure  of  the  saturated  liquid. 

3 

As  all  of  the  above  programs  assume  densities  in  g/cm  , temperatures 
in  kelvins  and  pressure  in  atmospheres,  a set  of  simple  routines  is  presented 
in  Appendix  III  for  conversions  back  and  forth  between  these  units  and  other 
commonly  used  systems  of  units. 
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Appendix  IV  contains  examples  of  typical  calculations  using  the 
routines  for  Appendices  I through  III. 

All  quantities  other  than  pressure,  density  and  temperature  are 
calculated  and  used  in  dimensionless  units.  Results  in  two  commonly  used 
systems  of  units  can  be  obtained  as  follows: 

To  obtain  heat  capacities  and  entropy  in  SI  units  (joules/g),  multiply 
the  dimensionless  quantities  by  0.48820;  or  in  Btu/lb.,  multiply  the 
dimensionless  quantity  by  0.210027. 

To  obtain  enthalpies,  internal  energies  or  heats  of  vaporization  in 
SI  units  (joules/g  K) , multiply  the  dimensionless  quantities  by  (0.48820  T) 
where  T is  in  K;  or  in  Btu/lb.  deg  F,  multiply  the  dimensionless  quantities 
by  (0.210027  T)  where  T is  in  deg  R. 

The  programs  as  presented  are  in  single  precision  except  for  the 
coefficients  of  the  generating  function  'Q',  which  are  in  double  precision. 
This  will  give  accuracy  to  within  the  limits  of  the  equation  when  used 
on  most  computers.  (When  used  on  computers  with  word  lengths  of  48  bits 
or  more,  even  these  coefficients  may  be  used  in  single  precision.)  If  more 
precision  is  needed,  the  following  can  be  made  double  precision  in  all 
routines:  T,  P,  D,  Q01,  Q02,  Q10,  Q20,  Oil.  The  greater  precision  may 
be  required  when  making  calculations  at  liquid  densities  or  in  the  vicinity 
of  the  critical  point,  particularly  when  obtaining  values  for  the 
thermodynamic  functions  requiring  higher  derivatives  of  the  generating 
function  Q,  such  as,  for  instance,  C^.  In  these  cases  the  contributions 
of  all  of  the  44  terms  are  of  importance,  as  there  is  much  cancellation 


of  contributions  from  individual  terms. 
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APPENDIX  I 

PACKAGE  OF  FORTRAN  PROGRAMS  TO  GENERATE  THE  THERMO- 
DYNAMI  C PROPERTIES  OF  AMMONIA  FOR  ENTIRE  REGION. 

SUBROUTINE  QQC Q ,L . K , T, D ) 

THIS  IS  THE  GENERATING  FUNCTION  FOR  Q AND  ITS  DERIVATIVES  USED  IN 
THE  CALCULATION  OF  THE  OTHER  FUNCTIONS.  Q IS  GENERATED  AS  A FCT 
OF  T AND  D.  AND  L IS  THE  ORDER  OF  THE  DERIVATIVE  OF  Q WITH  RESPECT 
TO  RECIPROCAL  TEMP.  AND  K IS  THE  ORDER  Oc  THE  DERIV.  W.R.T.  DENSITY. 
DOUBLE  PRECISION  A,  A 1 , A2 • D D.  C . TT.  QT 

DIMENSION  PC  5 » 5 ) , DDC 10>.Tr<7).AC44),IIC44).JJC44).AlC22)«A2C22) 
EQUIVALENCE  CA1C1).AC1)).(A2C1).AC23)) 

DATA  I 1/6*1 ,642,643.644,4*5. 4 *6, 4 *7 .4* 8.4*9/, N/4 4/ 

DATA  JJ/1.2.3.4.5.6.1.2.3.4.5.6.1 .2.3.4*5.6. 1. 2. 3. 4. 5. 6.1. 2.3.4. 1 * 

1 2.  3.  4.  1 . 2 . 3 . 4.  1 . 2 . 3 . 4 . 1 . 2 . 3 . 4 / 

DATA  P/2. . 2*. 3. . 4. . 5. . 2. .4. .6. .12. .20. .3. .6. .12. .24. *60. .4*. 12. * 

1 24. .48. .120. .5. .20. .60. .120. .240./ 

DATA  A 1/-6. 453022 30405D0. -13. 71 99  2 67705D0. -8.100  6 2031  57D0. -4 . 88009 
16421D0.— 12.0  28775626  0DO. 6.80634  593D0. 8.0800943676900. 1 4 . 3569 200 06D 
20. -45.0529767D0.-166. 1 8899  857D0, 3 7. 9 089  5023D0. -4  0. 730  2 083337DO. 

3 1 . 0329948807D0. 55. 84395581D0. 492 . 01 665081 77D0 . 1 737 . 8 36D0 . 

4 -30.87491 5263766D0. 7 1.48 3530 41 6D 0.-8. 948 264632D0, -16 9. 777744 1400. 

5 - 1 236. 532371 67D0. -781 2.  161 16  831 7D0/ 

DATA  A2/1 • 7795482691 4D0. -38.97461 096D0. -66 . 9 22858 82D0, -1.75394377532 
1 3 2D  0 . 208. 553371335D0.21348.9466D0.247.3417  46DO. 299. 98  3915550  0,4509 
2. 0 8057 879D0. -37980. 85D0,  -306.  55  78  854  300,24.  1 1 6 55 1 1 DO. -9323. 3 568D0. 
342724. 0985306D0, 1 61 . 791 0 0333746D0, -507. 478 0704 64 DO .3 1 39. 4703974D0 . 

4 -27458.71 063D0. -27.821 6888D0 , 298 . 81 291 73D0 .-2772. 597352D0. 

5 7668.92867800/ 

Q = 0 . 

IFCL+K)  12.14.18 
12  RETURN 

14  U*500./T 

C=U-1 .233349778D0 
IF(DABSCC) .LT.l.D— 8)  C=l.D-8 
TTf 1 >=1 . 

DO  15  1=2,7 

15  TTC  I 1 — TT C 1-1 »4C 
IFCD.LT.l.E-8)  D=l.E-8 
DDC 1 >=1 . 

DO  16  1=2.10 

16  DD( I I =DD C I — 1 >*D 
18  DO  200  M=1.N 

I=IIC  M) 

J=J J( M) 

IFCJ-L-1)  200,20.30 
20  I F ( L ) 200.24,28 
24  QT= 1 • 

GO  TO  50 
28  QT=P<L,L)/2. 

GO  TO  50 

30  IF  CL)  200.3  2,40 
32  QT=TTCJ) 

GO  TO  50 

40  QT=P<L. J-l )4TTC J-L) 

50  IFCK)  200.60.70 


c 


60  QT=QT*DD! I ) 

GO  TO  100 
70  R*l. 

DO  72  MM= 1 , K 
72  R=R*! I —MM ) 

QT=QT*R*DD! I-K) 

100  Q=Q+A!M)*QT 
200  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  FZ! T, FO. SO, CPO. UO ) 

C COMPUTES  THE  IDEAL  GAS  FUNCTIONS  AS  F CT  OF  T:  FREE  ENERGY! FO/RT ) 

C ENTROPY! SO/R) I CP!CPO/R>?  AND  INTERNAL  ENERGY! UO/RT » . 

DIMENSION  A ! 1 2 ) 

DATA  A/-3. 872727, .64*63724. 3. 2 238759. -. 0021 376925, .8689 0833E- 5. 

1 -• 240851 49E-7, .368931 75E- 10, -. 35034664E-1 3 . • 2 0563027E- 1 6 • 

2 -.665342E-20, .99392427E-24, 0 ./ 

U=T 

IFIT.LE.10.)  U=10. 

FO  = A < 1 l*ALOG! U)+A!2)/U+A! 3l+A!4l*U+A!5l*U*U+A!6)*U**3+A!7) *U**4 
1 +A!8)*U**5+A!9)*U**6*A! 1 0 )*J**7+A! 1 1)*U**8 

S0  = A!1  ) * ( A LOG!  U ) + 1 • ) + A!3>  + 2.*A!4)*U  + 3.*A!5)*U*U  + 4.*A!6)*l 

1**3  * 5.*A(7)*U**4  * 6 • * A < 8)*U**5  ♦ 7.*A!9l*U**6  ♦ 8. * A ( 1 0 ) *U** 7 
2+9. *A! 11 )*U**8 

CV0=A! 1 )/U  + 2. *A! 4) +6.*A! 5»  *U+1 2.*A! 6)*U*U*20 .*A! 7 ) *U* * 3 + 30. * A ( 8 >« 
1 u* *4  + 42. *A! 9)*U**5  ♦ 56 . *A! 1 0 ) *U**6  ♦ 72. *A! 1 1 ) *U**7 

CV0=-CV0*U 

SO=-SO 

UO  = SO+FO-1  • 

RETURN 

END 

C 

FUNCTION  PRES4D.T) 

C PRES  CALCULATES  THE  PRESSURE  IN  ATM  D SPHERES  USING  D.T.Q  AND  DQ/DlD). 
COMMON  /QQQQ/  QOO • QO 1 • QO 2 • Q 1 0 • Q20 , Q 1 1 
PRES  = 4.818*T*D*!l.+D*!Q00+D*Q0l  > ) 

RETURN 

END 

C 

FUNCTION  DPDD(D.T) 

C DPDD  CALCULATES  THE  DERIVATIVE  OF  THE  PRESSURE! IN  ATM)  WITH  RESPECT 
C TO  DENSITY  !IN  GM/CC). 

COMMON  /QQQQ/  Q 00 . QO 1 • QO 2 • Q 1 0 . Q20 , Q1 1 

DPDD  = 4.818*T*! 1 .+2.*D*Q00  ♦ D *D * ! 4 . *Q0 1 +D *Q0 2 ) > 

RETURN 

END 

C 

FUNCTION  ENTR! D. T. SO ) 

C CALCULATES  THE  ENTROPY  IN  DIMENSIONLESS  UNITS  !S/R).  REQUIRES 
C PRIOR  CALL  TO  FZ.  QOO  AND  Q10 

COMMON  /QQQQ/  Q 00 • QO 1 , QO 2 • Q 1 0 • Q20 • Q 1 1 

U=500 ./T 

TD=T*D 

IFITD.LT • 1 • E- 6)  TD=l.E-6 

ENTR=SIG  + D*U*Q10-D*Q00-ALOG<  4.818*TD ) 

RETURN 

END 


c 

FUNCTION  ENTH(D.T.ENER) 

C CALCULATES  THE  ENTHALPY  IN  D I ME MS  IONLESS  UNITS  (H/RT).  REQUIRES 
C PRIOR  CALL  TO  ENER,  FZ.QOO.QOl  AND  QIO. 

COMMON  /QQQQ/  QOO • QO 1 • QO 2 . Q 1 0 • Q20 . Q1  1 
ENTH=ENER  ♦ 1 • + D*QOO  «•  D*D*Q01 
RETURN 
ENO 
C 

FUNCTION  ENER ( D • T.  UO ) 

C CALCULATES  THE  INTERNAL  ENERGY  IN  DIMENSIONLESS  JNITS  IU/RT). 

C REQUIRES  PRIOR  CALL  TO  FZ,  QOO  AND  QIO 

COMMON  /QQQQ/  Q 00 • QO 1 . Q02 • Q 1 0. Q 20 • Q1 1 
U=500./T 

ENER=D*U*Q10  «■  UO 
RETURN 
END 
C 

FUNCTION  CV( D ,T  ,CPO I 

C CALCULATES  THE  HEAT  CAPACITY  AT  CONST  VOLUME  IN  DIMENSIONLESS 
C UNITS  (CV/R).  REQUIRES  PRIOR  CALLS  TO  FZ.  QOO  AND  Q20. 

COMMON  /QQQQ/  QOO • QO I • Q02 • Q 1 0 , Q20 • Q1 1 
U=500./T 

CV*CP0-1 .-D*U*U*Q20 
RETURN 
ENO 
C 

FUNCTION  CP ( D, T . CV ) 

C CALCULATES  THE  HEAT  CAPACITY  AT  CONST  PRESSURE  IN  DIMENSIONLESS 
C UNITS  CCP/R).  REQUIRES  PRIOR  CALLS  TO  FZ.  QOO.  Q01.Q02.  QIO.  Q20. 

C Q 1 1 AND  CV. 

COMMON  /QQQQ/  QOO • QO 1 • QO 2. Qt 0 • Q20 • Qt 1 
U=500./T 

CP^<  1 .«-D*(  Q00—U*Q1 0+D*  ( QO  1 -U*Q1 1 ) ) > **2/C  1 • ♦D*C  2 • *QOO*D*(  ♦ • *Q0  1*D* 
1 Q02 ) ) )+CV 
RETURN 
END 
C 

SUBROUTINE  OFI N0( DOUT , P.D. T. DPD) 

C THIS  SUBROUTINE.  BY  AN  ITERATIVE  PROCESS  USING  THE  - UNCTI ONS 
C PRES  AND  DPDD . W ILL  CALCULATE  THE  DENSITY  AS  A FCT  OF  A GIVEN 
C PRESSURE  AND  TEMPERATURE.  THE  REQUIRED  INPUTS  ARE  P.  T.  AND 
C DIWHICH  IS  AN  INITIAL  GUESS  FOR  THE  DENSITY! * THE  OUTPUTS  ARE 
C DOUT  CTHE  DENSITY  CONSISTANT  WITH  THE  INPUTS  » AND  T>  AND  DPD 
C (THE  VALUE  OF  DP/DD  AT  P AND  T).  THE  VALUES  OF  P ARE  IN  ATM. 

C T IN  K.  AND  D AND  DOUT  IN  G/CC.  NOTE  THAT  FOR  T LESS  THAN  TC. 

C AND  FOR  P AND  T NEAR  THE  COEXISTANCE  CURVE.  CARE  MUST  BE  USED  IN 

C THE  SELECTION  OF  THE  INITIA.  GUESS  FOR  D SO  THAT  SPURIOUS  VALUES 

C IN  THE  TWO  PHASE  REGION  ARE  NOT  RETURNED. 

COMMON  /QQQQ/  Q0.Ql.Q2.Ql0.Q20.Qll 

DD=D 

L*0 

9 L=L*1 

IFCDD.LE.O.)  D0=l.E-8 
I F < DD • GT . . 85 ) 0D=.85 
CALL  QQ(QO.O.O.T.DD) 

CALL  QQ( Ql.O.l.T*  DD ) 

CALL  QQ( Q2.0.2.T.DD) 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


DPDX=1 • 1*DPDDCDD.T) 

I F ( OP OX  . LT . • 1 ) DPDX=.l 
PP=PRESCDD,T) 

IFCABSC1, -PP/P ) . L T . 1 « E— 6 ) GO  TO  20 
X= { P— PP  J/DPDX 

I F( DABS( X ) • GT • • 1 ) X- XT • 1 /D A3S ( X ) 

DD=DD  4-X 

I F CL • LE . 20 ) GO  TO  9 
20  CONTINUE 
OOUT=DD 
RETURN 
END 

SUBROUTINE  PH ASEC P. T, DL • DG. Kl 
THIS  SUBROUTINE  W ILL  PERFORM  THE  FOLLOWING  FUNCTIONS: 

1.  IF  P AND  T ARE  SUPPLIED,  A VALUE  OF  K WILL  BE  RETURNED 

OF  —1,0, 1,2  OR  3 WITH  THE  FOLLOWING  MEANINGS!  K=-l  IS 
ERROR  RETURN  CP  OR  T LESS  THAN  0.1  K=0  POINT  <P,T)  IS  ON 

COEXISTANCE  CURVE  AND  APPROXIMATE  VALUES  FOR  DL  AND  DG  ARE 
CALCULATED,  K=1  INDICATES  T LESS  THAN  TC,  POINT  CP.T)  IS 
IN  VAPOR  PHASE.  APPROXIMATE  D=DG*DL  CALCULATED.  K*2  SIM- 
ILAR TO  K— 1 * BUT  POINT  ( P. T ) IS  IN  LIQUID  PHASE,  K*3 
INDICATES  T G.T,  TC,  AND  D=DG=DL  IS  APPROXIMATED. 

2.  IF  T LESS  THAN  TC  AND  P = 0,  PSCT)  WILL  BE  CALCULATED  ALONG 
WITH  APPROXIMATE  DG  AND  DL . IF  LESS  THAN  PC  ANO  T=0, 
TSCPI  WILL  BE  CALCULATED  A_  ONG  WITH  APPROXIMATE  DG  AND  DL. 

IFCP)  10,10.20 

10  IFCT.LT.195.  .OR,  T.GE.405.5)  GO  TO  90 
1 4 P=P S< T ) 

16  CALL  DSCT.DL.DG) 

K = 0 

RETURN 

20  IFCT)  22,22.30 
22  I F < P— 1 11,55)  24.24,90 
24  IFCP-.06)  90,26.26 
26  T=TSCP) 

GO  TO  16 

30  IFCT-405.5)  32.50,50 
32  PX=PSCT) 

34  CALL  DSCT.DL.DG) 

IFCP-PX)  36.14.40 
36  DG=DG*P/PX 
DL*DG 
K=1 

RETURN 
40  DG=DL 
K*2 

RETURN 

50  D=P*< .00021 07-5.19D-8*P-2.289D-7*T)-1.2823D-6*T 
D=D-*-P/T/4.81  8 
IFCD.LT.l.D-5)  D=. 00001 

DL  = D 
DG  = DL 
K = 3 

RETURN 
90  DL=0. 


! 


DG=0. 

K*-l 

return 


APPENDIX  II 

FORTRAN  PROGRAMS  FOR  FUNCTIONS  PERTAINING  TO  COEXISTANCE  LINE 
FUNCTION  PS(T) 

PS  CALCULATES  THE  VAPOR  PRESSURE  (IN  ATM)  AS  A FCT  OF  T IN  K. 

FOR  T GREATER  THAN  TC,  0.  IS  RETURNED, 

DIMENSION  A( 4 I/-7, 296510, 1 ,61 8053.-1 .956546 • -2.114118/ 

IF< T.GT. 405.4 ) GO  TO  10 

TH=T / 40  5,4 

TM= 1 ,-TH 

TS=DSQRT(TM) 

p*TM*TS 

Q=p*TM 

R — Q4Q 

F=  ( A ( 1 )*TM*A(  2)*P+A(3  )*Q+Af  4)  *R  )/TH 
PS=U1.85*EXP(F) 

RETURN 
10  PS=0. 

RETURN 

END 

FUNCTION  DPSDT(T) 

CALCULATES  THE  DERIVATIVE  OF  THE  VAPOR  PRESSURE  CURVE  WITH 
RESPECT  TO  T.  FOR  T GREATER  THAN  TC*  THE  VALUE  2.03  (THE  APPROX- 
IMATE VALUE  AT  TC)  IS  RETURNED. 

DIMENSION  A<  4) /-7.296510. 1 .618053.-1 .956545,-2.1 14118/ 

I F ( T • GT . 40  5.4)  GO  TO  10 

TH=T/405.4 

TM=1 .-TH 

TS=SQRT(TM ) 

p*TM*TS 

Q=P*TM 

R-Q4Q 

F— ( A ( 1 )*TM+A( 2)*P4A( 3)*Q^A(4)4R)/TH 
PS=11 1.85*EXP(F  > 

S=P*Q 

DPSDT  = ( — F — A ( 1 )— 1 . 5* A ( 2 ) *TS— 2 • 5* A ( 3 )*P-5 • *A  ( 4 ) * S ) 4PS/T 
RETURN 

10  DPSDT *2.03 
RETURN 
END 

FUNCTION  TSIP) 

CALCULATES  BY  AN  ITERATIVE  METHOD  USING  PS(T>  AND  GENERATING 
A STARTING  VALUE  BY  AN  APPROXIMATE  EQUATION.  THE  SATURATION  TEMPERATURE 
AS  A FUNCTION  OF  THE  PRESSURE  IN  ATM.  IF  P LESS  THAN  THE  TRIPLE  POINT 
VALUE  OR  GREATER  THAN  PC.  0 IS  RETURNED. 

Q=1 ./P 

t F( P . GT . .0  58  .AND.  P.LT.111.85)  GO  TO  4 
TS=0. 

RETURN 

4 CONTINUE 

T = 252.10703  ♦ P*(5. 6658636  ♦ P*  f -.12457077  *•  P*(. 00135987 
1 -P*. 00000531942  ))  ) Q* ( -1  5.  3 7967 84-Q* ( 1 . 52831  5-Q * • 0490 6 1 ) 

5 PX*PS(T) 

IF( ABSC l.-P/PX).LT.l.E-5)  GO  TO  10 


dt=cp-px)/dpsotit ) 

T=T+DT 
GO  TO  5 
10  T5=T 
RETURN 
END 

SUBROUTINE  DS<  T » DL, DG) 

C THIS  SUBROUTINE  WILL  CALCULATE  APPROXIMATE  COEXISTING  LIQUID 
C AND  VAPOR  DENSITIES  AT  THE  TEMPERATURE  SUPPLIED,  USEFUL  FOR 
C INITIAL  GUESSES  IN  THE  SUBROUTINE  0s  I ND,  IF  T GREATER  THAN  TC, 

C DL  = DG  = 0, 

DC= .235 

IFIT-405.4)  10*80,90 

10  DT=1 .-T/405.4 

IFIT.LT. 395. ) GO  TO  100 
T1*DT**.35 
T2*DT * * , 54 1 
T3=T2*T2 
T4=DT*2.9653 

DL=T4*T1*< 2. 117-1 .40974T2-. 83 802*T3 > 

DG=T4-T1*<2.1  174*1 ,1390*T2*.  57  253*  T3) 

DL*DC*(DL4-1  . ) 

DG— DC*(  DG-*- 1 . ) 

RETURN 
80  DL=OC 
DG=DC 
RETURN 
90  DL  = 0 • 

DG=0. 

RETURN 
100  Q=DT 

DL  = . 38  71  31  - • 00  09694  7/ Q 4- Q * I 1 . I 51  3 875  4-Q*  ( -1  • 4 9431  0 6 4-Q*  1.118332  5)) 
DG=.  082886  74-. 0 0095867/Q4-Q*  < - . 600  3953A4-Q*(  1 • 44  34  5 94-Q*  1 . 1678605) 
RETURN 
END 
C 

FUNCTION  H V ( T . DL, DG ) 

C THIS  FUNCTION  WILL  CALCULATE  THE  HEAT  OF  VAPORIZATION  IN  DIM- 
C ENS I ONLESS  UNITS  (L/RT)  AT  THE  INPUT  TEMPERATURE  T,  DL  AND  DG 
C ARE  REQUIRED  TO  BE  PREVIOUSLY  CALCULATED  AND  SUPPLIED  AS  INPUT 
U=500./T 

CALL  QQ ( QO . 0 , 0 • T , DL ) 

CALL  QQ(Q1 ,0,1,T,DL) 

CALL  QQ(Q2.1,0,T, DL ) 

CALL  QQ(R0.0.0,T,DG) 

CALL  QQ(R1«0.1«T, DG ) 

CALL  QQ1R2, 1 ,0,T,DG) 

HV=U*  ( DG#R2— DL*Q2  I ♦ DG*(  R04-DG*R1  ) - DL*  < QO 4-DL* Q1  ) 

RETURN 

END 

C 

FUNCTION  C S ( D • T , CPR ) 

C CALCULATES  THE  HEAT  CAPACITY  OF  THE  SATURATES  FLUID  IN 
C DIMENSIONLESS  UNITS  CCS/R).  INPUTS  ARE  O, T AND  THE  DIMENSIONLESS 
C CP/R  CALCULATED  PRIOR  TO  THE  CALL. 

COMMON  /QQQQ/  QOO.Q01,Q02,Q10,Q20.Q11 
TDPDT  = PRES(  D,  T ) — D*500  • *4.  81  8*D*  < Q1  04-D*Ql  1 ) 
CS*CPR-TOPDT/DPDD(D,T)/D/D  *D P SDT C T ) * . 1 0 1 325 
RETURN 


FUNCTION  TK(TF) 

TK=(TF+459.67)*5./9. 

RETURN 

END 

FUNCTION  TF(TK) 

TF  = 1 . 8TTK-459.67 

RETURN 

END 

FUNCTION  PATM(PSH) 
PATM=PSI A/14.696 
RETURN 
END 

FUNCTION  PSIA(PATM) 
PSI A=PATM*14.696 
RETURN 
END 

FUNCTION  P AT ( P 0 AR  ) 
PAT  =PBAR/ 1 .01325 
RETURN 
END 

FUNCTION  PBAR(PATM) 
PBAR=PATM *1.01325 
RETURN 
END 

FUNCTION  DGCC(QLRCF) 
DGCC=DLBCF/62 .428 
RETURN 
END 

FUNCTION  DLBCF(DGCC) 
DLBCF=DGCC*62 .428 
RETURN 
FND 


c 

C APPENDIX  IV 

C EXAMPLES  OF  USE  OF  PROGRAMS  IN  ABOVE  PACKAGES 

C 

C THE  FOLLOWING  SAMPLE  PROGRAM  WILL  CALCULATE  THE  DENSITY  AND  FIVE  THER 

C DYNAMIC  FUNCTIONS  AT  A POINT  <P,T)  SUPPLIED*  P SUPPLIED  IN  BAR,  T IN  DE 
DENSITY  IN  GM/CC,  CP.CV.S  IN  JOULES/GM  DEG  C.  AND  INTERNAL  ENERGY,  ENTh 
C IN  JOULES/GM. 

COMMON  /QQQQ/  QOO • QO 1 • QO 2 , Q 1 0 , Q20 • Q1 1 

1 READ  2,  PBAR.TC 

2 FORMAT  < 2F 10.3,F|0.6, 2F 10.3, 3F 10,5) 

P=PATM ( PBAR  I 

I F<  P.LE. 0,1  STOP 
T=TC+273.15 

CALL  PHASE  < P , T , DC.  , DG. K ) 

IF (K.EQ.-l ) GO  TO  1 
CALL  DF  I ND ( D » P * DL , T , D X I 
CALL  QQ(Q1 0, 1 ,0,T,DI 
CALL  QQ( Q20,2,0,T,D) 

CALL  QQ(Q11,1»1,T,D) 

CALL  FZ(T ,F0»  SO ,CV0,E0) 

CVV=CV(D.T.CVO > 

CPP=CP(D,T.CVV)*.4882 
CVV=CVV*.4882 
S=ENTR(D,T,SO )* .4882 
E=ENER  < O, T . EO ) 

H=ENTH<  D.T.E) *.4882*T 
E=E  * . 4882*T 

PRINT  2,PBAR.TC,D,H,E.S,CPP,CVV 

GO  TO  l 

END 

C 

C 

C A SAMPLE  PROGRAM  TO  CALCULATE  THE  SATURATION  TEMPERATURE 

C IN  DEG  F,  DENSITIES  IN  LB/CU  FT.  AND  HEAT  OF  VAPORIZATION  IN 

C BTU/LB  AT  A GIVEN  PRESSURE  IN  PSIA 

C 

REAO  1.  PSI 
1 FORMAT ( 5F 10.3) 

P = PAT  M ( PSI  ) 

T = TS<  P I 

CALL  DSCT, XDL.XDG ) 

CALL  DF I N D ( DL . P • X DL  » T , D X ) 

CALL  DFINDf DG,P. XDG.T.DY ) 

HH=HV(T,DL.DG»* .2100274T 
DL=DLBCF( DL > 

OG=DLBCF(DG> 

PRINT  1,  PSI ,TT,DL,DG.HH 

STOP 

END 
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